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Abstract —Measures of spike train synchrony have proven 
a valuable tool in both experimental and computational neu¬ 
roscience. Particularly useful are time-resolved methods such 
as the ISI- and the SPIKE-distance, which have already been 
applied in various bivariate and multivariate contexts. Recently, 
SPIKE-Synchronization was proposed as another time-resolved 
synchronization measure. It is based on Event-Synchronization 
and has a very intuitive interpretation. Here, we present a 
detailed analysis of the mathematical properties of these three 
synchronization measures. For example, we were able to obtain 
analytic expressions for the expectation values of the ISI-distance 
and SPIKE-Synchronization for Poisson spike trains. For the 
SPIKE-distance we present an empirical formula deduced from 
numerical evaluations. These expectation values are crucial for 
interpreting the synchronization of spike trains measured in 
experiments or numerical simulations, as they represent the point 
of reference for fully randomized spike trains. 

I. Introduction 

Understanding the details of information processing in the 
brain is one of the most challenging and exciting problems of 
our time. It has been widely established that the brain can be 
considered as an enormous network of spiking neurons, and 
that the spikes convey the information processed within this 
network ID. El- In the first studies, it was assumed that the 
information is encoded in the spike rates of the neurons, and 
several experiments confirmed this assumption 0. However, it 
became clear rather soon that in many, typically more complex 
situations the rate coding is not sufficient to represent the 
available information. For example, it was found that even 
single spike events can be responsible for the discrimination 
between different stimuli [01 • Hence, solely studying the spike 
rates is not sufficient for unraveling the neural code - exact 
spike timings need to be analyzed. 

In the last two decades various spike train distances, some 
inspired by existing mathematical distance measures, have 
been proposed 0, 0, 0, 0, QDJ, DU, DU, 03, DU¬ 

DS, DS- Their application to real neural data has led to a 
remarkable increase in the understanding of neural networks 
and neural coding na. 

An important step was the introduction of time-resolved 
synchrony measures that allow to analyze the time dependence 
of spike train similarities 0. lITOl . 114| . With this ability, it 
is now possible to investigate synchrony changes of pairs or 
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groups of neurons, for example induced by external or internal 
stimuli. Hence, such time-resolved synchronization profiles 
open new opportunities in spike train analysis. 

However, before proceeding with a detailed analysis of spike 
train similarity, first one has to obtain the spike events from the 
experimental measurements. This process is typically called 
spike detection and its difficulty is often underrated. Therefore, 
we give a brief overview of the existing techniques of spike 
and event detection in Section[EI| which represents the basis for 
the analysis of spike train synchrony discussed in Section [III] 
There, we present three time-resolved synchrony measures: 
the ISI-distance DID, the SPIKE-distance [Q3 as well as 
a new, time-resolved variant of Event-Synchronization 0, 
called SPIKE-Synchronization fTSl . The ISI- and SPIKE- 
distance, discussed several times in the past and being well 
formalized, still miss a comprehensive analysis of their math¬ 
ematical properties. Here, we will fill this gap and provide 
mathematical details of the ISI-distance in Section IIII-AI and 
the SPIKE-distance in Section Ull-B I Section ITlI-CI introduces 
SPIKE-Synchronization, including a detailed analysis of its 
mathematical properties. Finally, Section [IV] contains a brief 
summary and our conclusions. 

II. Detection Methods 

Experiments on brain activity, whether in vivo or in vitro , are 
typically done by placing electrodes in the region of interest. 
Most common techniques provide intracellular recordings of 
a membrane potential of a single isolated neuron (single-unit 
recordings) or measure the mean extracellular field potential 
generated by electrical activities of several (nearby) neurons 
(multi-unit recordings). The fundamental observables are the 
spikes, also called action potentials, emitted by the neurons 
and propagated to target nodes of the neural network via 
synapses. As the name suggests, these spike events are rapid, 
distinct maxima of the membrane potential. The prerequisite 
for any spike train analysis is the extraction of the spike times 
from the measured membrane potential. The set of spike times 
of a single neuron is then called the spike train, and is defined 
as: 

s = {fj, with ti<t i+1 , (1) 

where the f t represent the times of the spikes. Below we give 
a brief introduction of the current spike and event detection 
methods before proceeding to the measures of spike train 
synchrony. 
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Fig. 1: Panel A shows a spike detection scenario in which the 
dotted line is an adaptive threshold equal to three standard 
deviations of the signal (neuronal membrane potential). Panel 
B shows the detection of two different types of events - 
upward zero-crossings and local maxima above one standard 
deviation of the signal (EEG). 


A. Spike Detection and Sorting 

Spike detection means to detect discrete spiking events in 
continuous profiles of membrane or extracellular field poten¬ 
tial. Most commonly, the first step is to increase the signal to 
noise ratio by enhancing the spike waveform and reducing the 
noise. The simplest way of determining the spike locations is 
a threshold detector m (cf. Fig. QJV). Many spike detection 
techniques are based on an amplitude threshold with little or 
no preprocessing and either static or adaptive thresholds 11201 . 
ED, G2, (23). In many experimental situations where the 
extra-cellular field potential is measured, the profiles contain 
spike events from more than one neuron. Therefore, one not 
only has to identify the spikes, but also assign them into groups 
corresponding to the different neurons. 

In cases of intracellular recording or isolated cells however, 
one can assume that all spikes of one recording belong to a 
single neuron. There exist several techniques of spike detection 
for this situation, the most important ones are briefly described 
below. 

Static threshold detectors rely either on a single threshold 
for the detection of one edge (20) or on several thresholds to 
detect the slopes of the spike ED or other characteristics such 
as skewness, width, area under the peak, etc. Although simple 
thresholding is attractive for real-time implementations due to 
its computational simplicity, it is thought to be too sensitive 
to noise and often requires user input to set effective threshold 
levels (23). Overlapping spikes further reduce the efficacy of 
simple threshold detectors (24). 

Adaptive thresholds are dependent on the nature of the 
signal and are particularly useful in multiple channel record¬ 
ings, allowing optimization in several channels simultaneously. 
The threshold levels are often based on the signal’s standard 
deviation (25) or the standard deviation of the noise computed 


with respect to the median of the signal in which the spiking 
activity has the smallest contribution [261. 

Energy-based spike detectors first apply a nonlinear oper¬ 
ation to the signal before using static or adaptive threshold 
detection (27), (24). This nonlinear energy operator estimates 
the square of the instantaneous product of amplitude and 
frequency of a sufficiently sampled signal Il28ll and increases 
the separation of the spikes from the background noise. 

Often, one faces the more complicated situation where the 
measured signal possibly contains spike events from several 
neurons (multi-unit recordings). Then, additionally to spike 
detection also spike sorting must be performed. This involves 
isolating the neural signals and assigning each recorded wave¬ 
form to the neuron of origin (23) . Spike sorting techniques are 
based on the fact that action potentials recorded from the same 
cell tend to have a stereotypical spike shape determined by 
the cell’s morphology and biophysical properties, but also by 
its position relative to the recording electrode (29) . After the 
spike events have been detected using the techniques described 
above, spike sorting is usually done in two steps. 

The first step is feature extraction, where a number of 
significant features are obtained from the spike events, which 
allow to separate the different clusters afterwards. The most 
prominent method used for feature extraction is wavelet anal¬ 
ysis (26), (30) . which was shown to outperform the previously 
used principal component analysis (PCA) ED. 

The second step is clustering, where the spikes are assigned 
into different groups corresponding to different neurons. The 
clustering is based on the features obtained before. Typi¬ 
cal approaches include expectation maximization (32) , super- 
paramagnetic clustering (SPC) (26) or support vector ma¬ 
chines (SVN) ||33|| . 

It should be noted that despite the progress made in recent 
years, spike sorting remains a hard problem. Typical issues 
are the drifting of neurons, overlapping spikes and neurons 
with very similar spike waveforms. Increasing the number of 
observation electrodes or simultaneous intra- and extra-cellular 
recordings improve the results, but even then for experimental 
data it is usually impossible to verify if the spikes were 
assigned correctly. 

B. Event Detection 

Spike detection is the most prominent example for a 
transformation from continuous data to a discrete point pro¬ 
cess. This transformation becomes more complicated if the 
continuous time series do not contain any pronounced and 
stereotypical events such as spikes. In this case it is possible to 
use rather generic events such as zero crossings, local maxima 
and minima or any other particular feature characteristic of 
the signal (7). An example of an electroencephalogram (EEG) 
trace with two different types of events is shown in Fig. EJ 3 - 

The appropriate definition of the event is crucial since the 
information from the continuous time series is reduced to 
the time stamps of the events. It is not the synchronization 
between the time series as a whole that is evaluated, but rather 
the synchronization between the defined events only. Different 
choices of events can yield different results. 
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ISI: /(f) 

SPIKE: 5(f) 

SPIKE-Sync: C k 

Type of profile 

piecewise constant function 

piecewise linear function 

discrete function 

Relation between pro¬ 
file and overall measure 

Dj = T fl(t)dt 

D S = ±J S(t)dt 

SYNC= ±Z k C k , 

Dsync = 1 - SYNC 

Range 

0 < 7(t),D/ < 1 

0 < 5(f), < 1 

0 < C k , SYNC, D S ync < 1 

Value for identical 
spike trains si = S 2 

/(f) = 0, Dj = 0 

S(t) = 0, D s = 0 

Ck, SYNC = 1, D sync = 0 . 

Expectation value for 
Poisson spike trains 
with r = A 1 /A 2 

( Dl ) - (l+r-p + (l +r —1)2 

(D s ) w 1 - f e - h°s r ) 2 /s 

(SYNC) - r _4 r+2> 

(Dsync) = 1 - (SYNC) 


TABLE I: Overview of the properties of different spike train synchronization measures. 


III. Distance Measures 

To quantify the degree of synchronization of two spike trains 
si and S 2 , a distance measure I) is introduced to map the 
pair of spike trains into a positive number representing the 
differences between these spike trains. The normalization of 
D is arbitrary, but a sensible choice is to limit the distance to 
the interval [0,1], i.e.: 

D : {si, s 2 } i-A [0,1], (2) 


where D = 0 represents identical spike trains si = s 2 , while 
larger values denote a higher degree of dissimilarity. Many 
distance measures rely on one (or more) parameters that have 
to be chosen appropriately for the given spike trains and are 
usually connected to the typical time-scale of the spike events, 
for example the Victor-Purpura and van Rossum distances 0. 
0- Clearly, the existence of such a parameter introduces 
ambiguity in the analysis as its optimal value is a priori 
unknown, furthermore it is even unclear how to define such 
an “optimal” parameter value |34| . 

To overcome this problem, parameter-free distance mea¬ 
sures have been introduced, namely the ISI- and the SPIKE- 
distance as well as SPIKE-Synchronization described below. 
The ISI- and SPIKE-distance are based on time dependent 
distance profiles /(f), S(t), from which the overall spike train 
distance can be computed via simple integration, e.g.: 

As = ^ J S(t)dt. (3) 

The time T denotes the duration of the spike trains, i.e. the 
recording interval. Having such a profile 5(f) allows for a 
time resolved analysis of the spike train synchrony, which is 
important for example to detect synchronization triggered by 
external or internal events. 

So far, we have only introduced a bivariate distance, that is 
the distance between two spike trains si and s 2 . But usually 
one has to deal with multiple spike trains, e.g. from simul¬ 
taneous measurements of several neurons. Bivariate distance 
can be extended to multivariate distance by simply averaging 
the bivariate distances of all pairs of spike trains. Suppose we 
have N spike trains, then the averaged multivariate distance 
is computed as: 


N -1 N 


D a = 


_ 11 


N(N-l) ^ ^ 

n—1 m—‘ 


(4) 


n+1 


where D n m = D(s n ,s m ) is the distance between the spike 
trains s n and s m . As this average commutes with the time in¬ 
tegration, we can also readily introduce a multivariate distance 
profile: 

N -1 N 

Sa ^ = N(N — 1) ^ ^ Sn,m( ' 5 ' ) 

' ' n= 1 m=n -\-1 

where still D a = f S a (t)dt/T. 

This completes the framework of time-resolved spike 
train distances. In the following we will provide a de¬ 
tailed description of the 1ST and SPIKE-distance and discuss 
their properties. Additionally, the recently developed SPIKE- 
Synchronization is presented, another time-resolved measure 
based on the well-known Event-Synchronization (7). Despite 
also being time-resolved, SPIKE-Synchronization has some 
fundamental mathematical differences compared to the isl¬ 
and SPIKE-distance, as will be explained later. An overview 
of the mathematical properties of all three measures is shown 
in Table HI 


A. ISI-Distance 


The ISI-dissimilarity-profile /(f), introduced by Kreuz et 
al. QO), was the first parameter-free, instantaneous measure 
of spike train synchrony. It relies on the ratio between the 
concurrent interspike intervals of the two spike trains. 

1) Definition: Let {f^} be the spike times of the first spike 
train, then the interspike intervals are defined as: 


,(i) 


f (1) 

l i +1 




( 6 ) 


( 2 ) 

and similarly v- for the second spike train. To arrive at a 
time dependent profile, the sequences of interspike intervals 
are transformed into piecewise constant functions: 


v WM(t) = vU’<- 2) for ff ) ’ (2) <f<f« i (2) , (7) 

which for each interval [f», f*_|_i) takes the value Vi. Fig. [2] 
shows an example of two spike trains with the definition of 
j/( 1 ).( 2 )(i) at some time f. The ISI-profile is then given as the 
normalized absolute difference of the interspike intervals: 


m 


\v w (t) ~u (2) (f)| 

max{u( 1 )(t), ’ 


f e [o,T], 


( 8 ) 


The interval [0, T] again denotes the observation period of the 
two spike trains, i.e. 0 < f^ 1 ^ 2 ) < y However, there is a 
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potential ambiguity concerning the first and the last interspike 
interval if the first (last) spikes do not coincide with the start 
(end) time of the observation interval, e.g. 7 ^ 0. As an 

improvement to previous definitions Da, we here introduce 
an edge-correction for the ISI-distance to estimate the first and 
last interspike interval Da. Therefore, for the first interspike 
interval we use the maximum of the distance between the 
start of the observation interval and the first spike, and the 
first known interspike interval: v(t < t±) = max{fi,f 2 — fi}. 
A similar estimation is performed for the very last interspike 
interval. This makes the ISI-profile Eq. <© well defined for 
the whole observation interval t £ [0, T], As the interspike 
intervals are piecewise constant functions, also the ISI-profile 
is piecewise constant. The ISI-distance is then given as the 
integral over the whole profile, cf. Eq. #• 

2) Properties: From the normalization in Eq. (S), it is 
immediately obvious that the ISI-profile is bounded by 
0 < 1 (f) < 1 , and the value /(£) =0 is only obtained for iden¬ 
tical interspike intervals u^\t) = v^ 2 \t). Consequently, also 
the ISI-distance is bounded: 0 < Dj < 1. The ISI-distance is 
clearly symmetric: Di(si,s 2 ) = -D/(s 2 ,si) by construction. 
Furthermore, it can be shown that the ISI-distance also fulfills 
the triangle inequality: Dj(si, s 2 ) + Dj(s 2 , s 3 ) > Dj(si,s 3 ). 
A proof is sketched in Appendix B of (35). For identical spike 
trains it is obvious that the ISI-distance vanishes: Dj = 0. 
However, for two spike trains with constant and equal inter¬ 
spike intervals u^(t) = v ( 2 ) (f) = v, but with a global shift, 
the ISI-distance also evaluates to zero. Due to this degeneracy, 
the ISI-distance is only a pseudo-metric , but a full metric space 
can be recaptured by considering all degenerate spike trains 
with the same constant interspike interval, but overall time 
shifts as an equivalence class. 

Fig. [3] shows exemplarily a multivariate ISI-profile of 
TV = 50 spike trains. The change from the noise dominated 
spikes to increasingly synchronous events is captured quite 
well by the ISI-profile. In the beginning, we observe ISI- 
values very close to those expected for random Poisson spike 
trains, while the values drop significantly in the second half 
which indicates higher similarity. However, the ISI-profile is 
unable to detect the additional synchronous events within the 
random spike trains at the beginning. This is due to the fact that 
the ISI-profile only incorporates information about interspike 
intervals, and not about exact spike timings. 

For the ISI-distance it is possible to compute the expectation 
value (Dj) for two Poisson spike trains analytically. Note, that 
this expectation value depends only on the ratio of the rates of 
the two spike trains: r = Xi/X 2 ■ A straightforward calculation 
of this value is performed in Appendix [A] and gives: 


(Mr)) 


1 ( 1 
(1 + r) 2 + (l + r- 1 ) 2 ' 


(9) 


For two Poisson spike trains with equal rates, r = 1, we thus 
find (Dj) = 1/2, while in the limit where one spike train is 
much faster than the other one, r —> 0, 00 we have ( Dj) —> 1. 
This is visualized in Fig. [4] 


B. SPIKE-Distance 

The SPIKE-dissimilarity-profile S(t), first introduced 
in llT3ll and subsequently improved in H4l . provides a time- 
resolved distance measure that relies on the exact timings of 
spike events. 

1) Definition: The computation of S(t) is based on the 
four corner spikes surrounding the current time t : the pre¬ 
ceding spikes tp^’^ 2 \t) and the following spikes (t) 

of each spike train (cf. Fig. [2). The current interspike in¬ 
tervals can be expressed in terms of these corner spikes as 
j/( 1 )>( 2 )(f) = For each of the corner 

spikes, the distance to the closest spike of the other spike 
train is computed, e.g.: 

A t ( p\t) = min{|fp' ) -f- 2) |}, (10) 

l 


and similarly for Afp’ and ANote, that by definition 
A tp^p 2 \t) are piecewise constant functions. Fig. j^j shows a 
graphical representation of these definitions. These distances 
are then weighted by the distance of the corner spikes to the 
current time by the weighting factors: 

x <XU2) {t)=t _ t DH2) and x (iU2) {t)=t (i),w_ tA id 

with + Xp' l ' l ' 2 \t) = j/( 1 )’C 2 )(f), as seen from Fig. [ 2 ] 

The weighted distance for the first spike train then reads: 


Si{t) 


Atp\t)xp\t) + Atp\t)xp' > 
!/(!)(*) 


( 12 ) 


and similarly S 2 (t) is defined for the second spike train. As 
xJp' l p 2> (t) are linear in £, the resulting functions are 

piecewise linear, with possible jumps at the interval edges 
£.(t).( 2 ). Finally, these local distances are then weighted by 
the local interspike intervals and with a proper normalization 
we arrive at the definition of the SPIKE-profile: 


Sit ) 


Si(£)z/ 2) (f) + S 2 (t)u ( ' 1 \t) 
l(yD)(t) + vM(t)) 2 


(13) 


which again is a piecewise linear function as S\ t2 are piece- 
wise linear while the other terms are piecewise constant. 

As for the ISI-profile above, we face a potential ambiguity 
for the first and the last interval. Previous definitions of 
the SPIKE-distance introduced auxiliary spikes at the edges, 
which lead to spurious synchrony. An improved treatment of 
the edges, similar to the edge correction for the ISI-distance 
above, is presented in fl8l . 

2) Properties: The normalization in Eq. m again en¬ 
sures that the SPIKE-profile is bound to 0 < S(t) < 1. 
Hence, the same bounds hold for the SPIKE-distance 
0 < Dg = f S(t)dt < 1. Furthermore, Ds(si,s 2 ) = 0 only 
if Si = S 2 , and the SPIKE-distance is also symmetric 
Ds(si,s 2 ) = Ds(s 2 , s\). However, it is not a metric as it 
is possible to construct spike trains that violate the triangular 
inequality: D s (s 1 ,s 2 ) +D s (s 2 ,s 3 ) ^ L> s (si,s 3 ). 

Fig.[3]shows an exemplary SPIKE-profile for the same TV = 
50 spike trains as before. As for the ISI-profile, the overall 
behavior of noise dominated spike times in the beginning and 
the increasing synchronization afterwards is well captured by 
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Fig. 2: Local definitions of interspike intervals and time 
differences required for the calculation of the ISI- and the 
SPIKE-profile. 



Time 

Fig. 3: Multivariate ISI, SPIKE and SPIKE-Synchronization 
profiles for M = 50 spike trains (shown in top panel). 
The spike trains are generated artificially such that the first 
half consists of a noisy background with a few synchronous 
events with increasing jitter, while the second half contains 
increasingly synchronized events. The dashed lines represent 
the respective expectation values for Poisson spike trains. 


the SPIKE-profile. Additionally, the SPIKE-distance, relying 
on exact spike timings, is able to detect the synchronization 
events in the beginning as seen from the clear minima in the 
SPIKE-profile. 

Unfortunately, the complicated definition of the SPIKE- 
distance makes an analytic computation of the expectation 
value ( Dg) for two Poisson spike trains intractable. However, 


it is clear that also the SPIKE-distance should only depend on 
the rate ratio r = A 1 /A 2 of the two Poisson spike trains and 
approach the value Dg = 1/2 in the limit r —► 0, 00 . This 
is seen in Fig. [4] Similarly to the ISI-distance and SPIKE- 
Synchronization below, the SPIKE-distance exhibits a clear 
minimum for spike trains with equal rates Ai = A 2 , i.e. r = 1. 
As we are unable to obtain an exact analytic result for ( Dg ) at 
this point, we provide an empirical approximation. Therefore, 
we use the following function that already incorporates the 
properties mentioned above (limr^o.oo (Ds) = 0.5, minimum 
at r = 1): 

(D s ) = * - ae- (logr)2 / (2/32) . (14) 

From visual inspection, we find that with a = 0.2, 0 = 2, 
Eq. ( [T4| provides an excellent approximation of the average 
SPIKE-distance for two Poisson spike trains as shown in 

Kg- E 


C. SPIKE-Synchronization 

SPIKE-Synchronization can be understood as an in¬ 
stantaneous coincidence detector. It was recently pro¬ 
posed by Kreuz et al. and is derived from Event- 

Synchronization (7). In contrast to the ISI- and SPIKE- 
distance, SPIKE-Synchronization quantifies similarity instead 
of difference, but a distance measure can be constructed in a 
straightforward way as shown below. 

1) Definition: For the SPIKE-Synchronization profile, a 
coincidence indicator is defined for every spike of 

the two spike trains This coincidence indicator can 

have two possible values: Cj = 1 if the spike at t, is part 
of a coincidence, and 6/ = 0 if not. Similar as for Event- 
Synchronization, this coincidence indicator is given by: 



1 if min^df. 1 ^ 
0 otherwise. 



(15) 


That means, a coincidence is found for the spike i of the first 
spike train if the distance to the closest spike j of the second 
spike train is smaller than the coincidence window r// . The 
coincidence window is defined adaptively according to the 
local firing rate: 


r (L2) 

1 ij 


,(!) „( 1 ) „( 2 ) ( 2 ) 


V' > \ 

v i- 1J> 


(16) 


with i/LA 2 ) being the interspike intervals as given in Eq. (|6b. 

(2Y"T 

The coincidence indicator for the second spike train C) is 
computed in the same way as in Eq. © but with exchanged 
indices (1) -f-)- (2). 

The SPIKE-Synchronization profile is then obtained 
by merging the coincidence indices of the two spike 
trains {C/} = {C'| 1) } U {C- 2) } as well as the spike times: 
{t' k } = {fp*} U {fj 2 *}, which results in a discrete function 
defining a SPIKE-Synchronization value for each spike time 
of the two spike trains: t' k 1 —>■ Cf . For visual purposes, it 
might be reasonable to connect the individual points of this 
discrete profile, but opposed to the ISI- or SPIKE-profile 
above, SPIKE-Synchronization is only defined exactly for 
the spike times - any connecting lines do not represent real 
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Fig. 4: Numerical results for the ISI-distance (circles), 
SPIKE-distance (triangles) and SPIKE-Synchronization- 
distance (squares) of two Poisson spike trains with a total 
of Mr: 20000 spikes in dependence of the rate ratio 
r = A 1 /A 2 . The black dashed lines represent the analytical 
result for the ISI-distance Eq. the empirical curve for 
the SPIKE-distance Eq. ( [T4| and the analytical result for the 
SPIKE-Synchronization distance Eq. (J22|i. 



intermediate values but are merely a guide to the eye. Note 
that Ck = 1 for all k means that every spike is part of a 
coincidence, while Ck = 0 means no synchronous spikes 
appeared. 

Being a discrete function, integrals are not defined for the 
SPIKE-Synchronization profile. However, an overall SPIKE- 
Synchronization value SYNC can be obtained very naturally 
by summation: 

1 M C 

SYNC = — V C k = —, (17) 

M ^ M 

k=1 

where M is the total number of spikes in the merged spike 
train {t' k } and C denotes the total number of coincident 
spikes. As seen from Eq. 0 . this value has a very intuitive 
interpretation as it simply represents the fraction of coinciding 
spikes. The coincidence factor introduced in li36l and further 
analyzed in lfl2l has a similar interpretation. However, it is not 
time resolved and its coincidences are defined with a constant 
coincidence window r = const. Note, that SYNC quantifies 
similarity, but a distance measure can be trivially obtained: 


Aync = 1 - SYNC. (18) 

The discrete nature of the SPIKE-Synchronization profile 
also requires a different definition for the multivariate case, 
as discrete functions can not be added in a straightforward 
manner as done for piecewise constant or piecewise linear 
functions in Eq. <0- For the multivariate profile of N spike- 
trains, we first define a generalized bivariate coincidence 
indicator C k "'’" ' for all pairs of spike trains: 

c (n,m) J 1 if minjfl t\ -t) J |) < T-y (19) 

1 0 otherwise, 


where rjj !1 is the same as defined in Eq. |l 6 j ), but for arbi¬ 
trary spike trains n and m. From these bivariate coincidence 
indicators we then compute the average coincidence counter 
for each spike in every spike train: 

c\ n) = £ <?("-"*). ( 20 ) 

m^n 


Finally, again all the averaged coincidence counters are merged 
{Ck} = U n {Ci™Y Together with the merged spike times 
this results in the multivariate SPIKE-Synchronization profile. 
Note that in contrast to the bivariate case, the multivariate 
profile can obtain values different from just zero and one. 
Namely, from Eq. (20 1 we find that C k = p/(N — 1) 
with p = 0... N — 1. Furthermore, note that for N = 2 
the multivariate definition in Eq. ( |20| becomes equivalent to 
the bivariate definition of Eq. (| 1 5|). Similar to the bivariate 
case, the overall multivariate SPIKE-Synchronization value is 
calculated as the ratio of coincident spikes: 


SYNC“ = 


A? 

M a 


C a 

W a 


A y nc = 1 - SYNC a , (21) 


where M a is the overall number of all spikes in the N spike 
trains. 

2) Properties: In contrast to the 1ST and SPIKE-profiles 
that quantify dissimilarity, the SPIKE-Synchronization pro¬ 
file is a similarity measure. Furthermore, the SPIKE- 
Synchronization profiles are defined on discrete points only 
(the spike times) and are discrete valued as well. While 
a bivariate profile can only take two values: zero and one 
representing the absence or presence of a coincidence, the 
multivariate profiles can exhibit N values in [0,1]. It is 
immediately clear that the SPIKE-Synchronization distance is 
not a metric, as there exist spike trains with -D S ync(si, S 2 ) = 0 
even for si / S 2 - Furthermore, also transitivity is violated 
as one easily finds examples where T , sync (s 1 , s^) = 0 and 
Tfsync ip‘2- £ 3 ) — 0, but D sync (si,s 3 ) > 0 . 

An exemplary multivariate SPIKE-Synchronization profile 
is shown in Fig. [3] As seen there, SPIKE-Synchronization can 
clearly differentiate between the noise dominated part at the 
beginning and the synchronous spikes at the end. However, at 
the transition in between, large fluctuations appear due to the 
discrete nature of the coincidence measure. Nevertheless, the 
synchronization events at the beginning are clearly captured 
as distinguished peaks. 

An analytical calculation of the SPIKE-Synchronization 
distance (D s ync ) of Poisson spike trains is performed in 
Appendix [B] and gives: 

(A y nc) = l- r ,_ 1+ 1 2 + r , (22) 


where r = A 1 /A 2 is again the ratio of the rates. 


IV. Summary and Conclusions 

In this paper, we first provide a comprehensive, concise 
introduction to the different techniques of spike detection. This 
is the basis for any kind of spike train analysis which can 
help to understand the fundamental processes in the brain. 
One approach is to analyze spike train synchrony and for this 
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purpose we here review three time-resolved measures: the isl¬ 
and the SPIKE-distance, which have been known for several 
years, and the recently proposed SPIKE-Synchronization. All 
these methods provide a way to quantify the similarity of 
spike trains in a local manner, but can also be reduced to 
an overall distance measure: Dj, D$ and D sync . Although 
originally defined only for bivariate profiles of two spike trains, 
they can be generalized to multivariate situations and measure 
the combined synchrony of an ensemble of spike trains. 

Furthermore, we study the mathematical properties of these 
three measures in order to provide the necessary information 
for their successful application to experimental or numerical 
data. Specifically, we are able to obtain the expectation values 
of the overall distances for two Poisson spike trains of arbi¬ 
trary rates. This gives another, crucial point of reference for 
interpreting experimental and numerical results. For the ISI- 
distance as well as the SPIKE-Synchronization, we can cal¬ 
culate the expectation values exactly, cf. Eq. and Eq. ( [22| . 
For the SPIKE-distance, the analytic result is intractable at this 
point and instead we present an empirical estimate in Eq. © 
that shows excellent agreement with numerical results, cf. 
Fig- El Note that for all three methods, the expectation value 
for two Poisson spike trains only depends on the ratio of the 
rates r = A 1 /A 2 . This is an immediate consequence of the 
invariance of the measures under global rescaling of the time. 

This detailed mathematical analysis of the spike synchro¬ 
nization measures further improves the usability of these meth¬ 
ods for examining and interpreting experimental and numerical 
results. 

We finally note that all three methods are part of the 
Matlab based graphical user interface SPIKyQ [HI, 1371 . 
Furthermore, with the PySpike library^] there also exists a 
Python implementation for spike train similarity analysis that 
provides the three methods discussed here. 
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Appendix 

A. ISI-Distance for Poisson Spike Trains 

We are interested in the average ISI-Distance (D/) of two Poisson 
spike trains with rates Ai and A 2 . We start by noting that the average 

1 http://www.fi.isc.cnr.it/users/thomas.kreuz/Source-Code/SPIKY.html 

2 http://www.pyspike.de 


ISI-distance is identical to the average profile: (Di) = (I[t)). 
Furthermore, for Poisson spike trains with rate A the probability 
density for the interspike interval v is P{v) = Xe~ Xv . However, if we 
consider the interspike intervals as a piecewise constant function v(t), 
cf. Eq. |7J, the probability to be in an interval of length v{t) at time t 
is the probability of the interspike interval multiplied by its length: 


P(u{t)) = A ~u(t)e 


—Xu(t) 


(23) 


where the extra factor A ensures normalization. For the average 
profile value we consider the interspike intervals as independent 
random variables vi t2 and perform the integration over the probability 
density: 


(I) = 


poo po 

/ dz/1 / 

Jo Jo 

POO P o 

/ di/1 / 

Jo Jo 


dv 2 I(ui, v 2 )P{v 1 )P{v 2 ) 


dv 2 


\Vi - V 2 \ 


\ 2 —Aizyi \ 2 —\ouo 

-Aiz/ie A 2 ^ 2 e 


maxjTi, i/ 2 } 

The absolute value and the maximum can be resolved by considering 
ui > v 2 and < o 2 separately and splitting the integral: 

r vi 


(I)=xl\l 


POO Pl> 

/ dvl / 

Jo Jo 


6 .U 2 (yi — V 2 )V 2 Z 


— — A 2 IY 2 


Xi 


POO PC 

XlXl / din / 

JO J v 1 


dv 2 {v 2 — vi)vie 


— A11/1 — A 2 ^2 


Ai 


The first integral evaluates to: 

rr _ -1 rj Al ^ Aj_ 

A 2 A 2 (Ai T A 2 ) (Ai + A 2) 2 ’ 
which can be expressed in terms of the rate ratio r = A 1 /A 2 : 


Xi = 1 - 2r + 


2r 2 


+ 


1 + r (1 + r) 2 (1 + r) 


2 ' 


(24) 


(25) 


For the second integral, after a change of the order of integration we 
obtain the same result as above but with 4-1 v 2 interchanged, hence 
with the inverse ratio r~ = A 2 /A 1 . Putting these results together we 
arrive at the average ISI-distance for two Poisson spike trains with a 
rate ratio r = A 1 /A 2 (see Fig. [ 4 J: 

1 1 


(Dj) = (1(f)) =T 1 +T 2 = 


+ 


(1 + r) 2 (1 + )— J ) 2 


(26) 


B. SPIKE-Synchronization Distance for Poisson Spike Trains 

Here, we calculate the expectation value of the SPIKE- 
Synchronization (SYNC) for two Poisson spike trains with rates Ai 
and A 2 . We start from the coincidence counter C \ 1 ] of the first spike 
train Eq. 1 15 j Note that it is sufficient to consider only C \' ^, as 
there we already account for all possible coincidences. The overall 
SPIKE-Synchronization can be calculated as two times the sum of 
this coincidence counter. For the average one thus finds: 


Mi 


(SYNC) = ( 


Mi + M 2 


E^> = 


2Ai 


Ai + A 2 


(Cf) 


(27) 


where Mi ,2 are the number of spikes in the spike trains and in the 
last step we used that (M) . 2 ) = Ai _2 for some spike train length T. 

To compute the remaining average, we express Cf in terms of 
the following three independent random variables: 


t = miriOf \ i/fi, tzfi}, P(t) = Ae At 


1 = 


P(y) = Alt/e"* 2 " 
P(x) = 1. 


(28) 


3 Choosing the first spike train is arbitrary, equivalently the second spike 
train could be analyzed with the same result. 

















The crucial point is that vf is removed from the minimum in t and 
considered separately. For Poisson spike trains the above variables 
can be interpreted as follows: for each spike f^ 11 we randomly choose 
the minimum of the three surrounding interspike intervals from an 
exponential distribution with rate A = 2 Ai + A2. Furthermore, the 
interspike interval v is chosen at random, where we have to take 
into account that the probability of finding some interval v around 
time is proportional to v. Thus we get P(v) ~ ve~ X2V . Finally, 
x € [0,1] is chosen uniformly to determine the distance between t 
and if. The coincidence counter is then expressed in terms of these 
random variables as: 

( 1 if xv < f/2 

Cf) = J ° r t 1 "*)^/ 2 (29 ) 

or v < t 
0 otherwise. 


The average is obtained from integrating over all random variables 
(Cf>) = f C^dfdi'do;. Substituting Cf, the integration over x 
can be performed immediately. The remaining integral yields: 


poo 

<cf >) = / dfP(f) 

Jo 

A 2 

“ 2 (Ai + A2) 

Consequently, using Eq. ( |27| i we find for the average overall SPIKE- 
Synchronization of two Poisson spike trains: 


[ d vP{v)+ [ d v-P{v) 

Jo Jf v 


(SPIKE) 


Ai A2 

(Ai + A2) 2 


1 

r + 1 — 1 + 2 ’ 


(31) 


where again the rate ratio r = A1/A2 is introduced and the result is 
symmetric in r and r -1 . Finally, the SPIKE-Synchronization distance 
amounts to: 

( Ay „c> = 1 - (SYNC) = 1 - r + r l 1+2 ■ (32) 

This result is numerically validated in Fig. [4] 
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